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Abstract 

Spiral surface growth is well understood in the limit where the step motion 
is controlled by the local supersaturation of adatoms near the spiral ridge. 
In epitaxial thin-film growth, however, spirals can form in a step-flow regime 
where desorption of adatoms is negligible and the ridge dynamics is governed 
by the non-local diffusion field of adatoms on the whole surface. We inves- 
tigate this limit numerically using a phase-field formulation of the Burton- 
Cabrera-Frank model, as well as analytically. Quantitative predictions, which 
differ strikingly from those of the local limit, are made for the selected step 
spacing as a function of the deposition flux, as well as for the dependence of 
the relaxation time to steady-state growth on the screw dislocation density. 
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Spiral surface growth is one of the most widespread growth mechanisms for crystals 
with atomically flat surfaces. Such crystals grow by the incorporation of new atoms at 
monoatomic steps. If steps are pinned at a screw dislocation, they wind around the dis- 
location and form growth spirals. The step dynamics, and the final steady-state spacing / 
between successive steps (or equivalently the surface slope a/ 1, where a is the lattice param- 
eter) is determined by the interplay of surface diffusion, the attachment kinetics of atoms 
at the steps, and the step line tension. Recently, there has been renewed interest in spi- 
ral surface growth following the observations of spiral ridges in sputtered high-temperature 
superconducting thin films |TJ and in certain semiconductor materials grown by molecular 
beam epitaxy (MBE) 0. 

In the classical Burton-Cabrera- Frank (BCF) model of surface growth atoms are 
first adsorbed to the crystalline surface ( "adatoms" ) and then diffuse along the surface until 
they are either incorporated into the crystal at a step, or desorb from the surface with a 
probability 1/t s per unit time. Therefore, two different growth regimes can be distinguished 
depending on the ratio of / and the diffusion length x s = y/Dr s , where D is the surface 
diffusion constant. When desorption is fast (x s <C I), only adatoms which are deposited 
near a step are incorporated, and the dynamics of the steps is local, that is the velocity 
of a step is completely determined by the local supersaturation and the step curvature. 
This is the well understood regime described by the BCF theory of spiral growth ||[|]. In 
many practical applications such as MBE, however, spirals can form in a step-flow regime 
at temperatures where desorption is negligible. Then, all deposited atoms reach a step, 
and successive turns of the spiral are strongly coupled via diffusion. The step dynamics 
becomes a highly non-local free boundary problem. The case of steady-state growth has 
been investigated by approximate theories |||| and by the boundary integral (potential 
theory) method [0], but no dynamical solution of the original equations has been performed 
to validate these results. Another important difference between the two regimes is their 
approach to steady-state growth: in the local limit, the spiral finds its final step spacing 
essentially after a single rotation. In contrast, without desorption one would expect a slower 
relaxation to steady-state due to the global redistribution of adatoms. In particular, this 
relaxation may depend on the density of screw dislocations. 

The main goal of the present letter is to use a phase-field approach to study the dynamics 
of spiral ridge formation in the BCF model. This approach has recently been used to effi- 
ciently solve a similar free boundary problem for dendritic growth || , and the mathematical 
results of this study are exploited here. One distinguishing feature of our approach is that it 
makes it possible to investigate the full crossover from the local to the desorpt ion-free limit, 
whereas previous works |],|nj have assumed a constant effective supersaturation, as appro- 
priate in the local limit. The phase field method in this context can also be interpreted as a 
direct continuum analog to microscopic growth models studied by Monte Carlo techniques 
[DJ . We make quantitative predictions for the selected step spacing as a function of the 
deposition flux and for the time to approach steady-state growth, under the assumption that 
the elastic interaction between steps can be neglected. We focus mainly on the situation 
where adatoms feel the same barrier for attaching at ascending or descending steps, i.e. no 
Ehrlich-Schwoebel (ES) barrier [T2| , in which case attachment can be described by a single 
sticking coefficient k. 
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We write the BCF equations in terms of the dimensionless diffusion field u = Q(c — c® q ), 
where c is the adatom concentration, Q is the atomic area of solid, and c° eq is the equilibrium 
concentration at a straight step. The basic equations have the form 

— = DV\--+F, (1) 




(2) 
(3) 

where v n is the step normal velocity, (du/dn)± is the normal concentration gradient on the 
lower (+) and upper (— ) side of the step, k is the local step curvature, d = f2 2 c° 7/(fcgT) 
where 7 is the step stiffness. The effective deposition frequency F is related to the actual 
deposition flux per atomic area F d by F = F d — c eq Q/r s . 

We simulate Eqs. 1-3 by reformulating them in terms of a phase-field model similar to 



the one used previously by Liu and Metiu [I3| for a one-dimensional step train. The basic 
equations of our model are 

dib 5H , , 

r ^ = -JI = ^ 2 VV + sin(7r[^-V s ]) 

+ A«[1 + COS(7T[^-'0J)], (4) 

where H is the free energy functional depending on the fields ib and u, ib/2 represents the 
surface height in units of a, ips/^ is the height of the initial substrate surface, u is the 
concentration field defined above, W is the width of the step, t$ is the characteristic time 
of attachment of adatoms at the steps, which is typically much smaller than r s , and A is a 
dimensionless coupling constant. If we replace the coupling between the two fields in Eq. |4] by 
a constant supersaturation, we obtain the continuum limit of the solid-on-solid model. The 



main difference with the model of Ref. fl3| is that the term [1 + cos(7r [i/j — ip s ])] is introduced 
in Eq. |] to keep the minima of H at fixed values (ip — ip s = 2n + 1), independently of the 
adatom concentration. A screw dislocation is introduced at the origin by choosing irip s equal 
to the polar angle in the x-y plane, or nips = atan(y/x), which corresponds to shifting the 
energy minima up by one atomic spacing after one complete counter-clockwise rotation. 

We now use the recent asymptotic analysis of Karma and Rappel M to relate the equa- 
tions of the phase-field model to the sharp-interface equations of the BCF model. Since the 
analysis of Ref. applies directly to Eqs. [| and|5|, it need not be repeated here. The results 
of interest are that in the thin-interface limit W/x s — > 0, the phase-field equations reduce 
to the free boundary problem defined by Eqs. |TJ and |^, with Eq. 3 replaced by 

u = d K + v n /k (6) 

where d = aiW/X and 1/k = a\{j^j (XW) — a 2 W/D). The kinetic coefficient k can then be 
made effectively infinite (instantaneous attachment), thereby recovering the condition ([3D 
of local equilibrium, by choosing the coupling constant A = t^D/^W 2 ). The numerical 
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FIG. 1. Step spacing I normalized by cIq as a function of F for finite and vanishing desorption. 
The standard BCF theory in the local limit is shown as a dashed line for comparison. 

constants a± and are determined by the form of the energy function in Eq. [| and can 
be evaluated using Eqs. 51, 54-56, and 58-59 in Ref. J5|, together with the one-dimensional 
stationary profile solution of Eq. | with u = for an isolated step with ^(=1=00) = ±1, 

1 / \ -, 4 / \Fkx\ 

ipo{x) = 1 atan exp ■ (7) 

71 \ W J 

The resulting numerical values are a\ = 0.718348 and 0,2 = 0.510442. Eqs. [| and [5] were 
simulated on a square lattice of edge length L = NAx with zero flux boundary conditions and 
iV varying between 100 and 600. As in Ref. ||, Eq. |] was integrated using an explicit Euler 
scheme, Eq. |5| using a Crank-Nicholson implicit scheme. For the simulations, we measured 
lengths and times in units of the phase-field parameters, i. e. W = = 1. In these units, we 
chose D = 10, which yields A = 19.591, d = 0.0366, Ax = 0.5, and At = 0.025. The final 
results, stated in dimensionless ratios of physically well-defined quantities, do not depend 
on the choice of W and as long as W <C x s , W <C /, and Tf 1/F. 

Simulations were started from a straight ridge along the x-axis with one end pinned to 
the screw dislocation (?/> = ip s ). After a long transient, the ridge reaches a steady-state spiral 
shape with a constant angular rotation frequency uj, and a step spacing /. For our plots, 
we defined / to be the distance between the first two successive steps from the center. This 
value is only a few percent smaller than the asymptotic step spacing far from the core. A 
plot of I as a function of F is shown in Fig. for x s /d = 136.4, and x s = 00 (no desorption). 

Let us now compare our results to the standard BCF theory of spiral growth that is 
based on the expression for the normal step velocity 

V n = foo (1 - r c K) (8) 

where v ^ is the velocity of a straight step and r c is the critical radius for island nucleation. 
Expressions for u and / are then obtained by looking for shape preserving solutions of Eq. 
p| in a frame rotating at angular frequency lu. In the local limit, I 3> x s , the solution is M: 
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I = 

it) = 



(9) 
(10) 



with = 2Fx s and r c = d /(Fr s ); Fr s is the supersaturation far away from steps. The 
corresponding curve for I is plotted as a dashed line in Fig. [I]. For an arbitrary ratio l/x s , 
Eq. [8] is no longer exact because the diffusion fields of the steps overlap. Cabrera and 
Coleman (CC) proposed || that a rough estimate of I and lo can be obtained by assuming 
that Eqs. PHTOl continue to hold with 

v O0 = 2Fx s t&vh{l/2x s ), (11) 
r c = d o /u(0) = d / [Ft s (1 - l/I (l/x a ))] , (12) 

where Eq. |TT] is the exact expression for the velocity of an infinite one-dimensional step 
train of spacing I, and «(0) is the supersaturation at the center of a circular terrace of 
radius I, obtained by solving Eq. [I] subject to the boundary condition u(l) = 0. This yields 
u(r) = Ft s [1 — Io(r /x s )/ Iq(1/x s )], where Iq is the zero-th order modified Bessel function. 
One interesting consequence of the CC estimate is that in the limit x s /l — > oo 

I = A(d D/F) 1/3 (13) 
u = 2nF (14) 

where A = 76 1/3 « 4.236. The expression for u) is exact and follows from global mass 
conservation: since for I <C x s all adatoms reach a step, the spiral rotation frequency must 
just be F in steady-state. The same scaling relations, but with a different prefactor A have 
been obtained by the boundary integral method 0. The scaling for I becomes exact in the 
limit do/ 1 — > 0, which is practically always satisfied. This can be understood from a simple 
dimensional analysis. In the limit where r s — > oo, all parameters can be removed from the 
BCF equations by making the variable transformations 

t' = tF, x' = x{F/d D) 1/3 , v! = u{D/d 2 F) 1/3 , (15) 

if one neglects dtu in Eq. 1 which is of relative magnitude do/ 1- Our numerical results 
confirm this scaling and yield a value of A = 4.626. They also show that the cross-over from 
I ~ F^ 1 to / ~ _F~V 3 is very slow. 

Next, let us examine how the dynamical aspects of the spiral ridge formation depend on 
the system size L, which on a real surface is roughly given by the mean distance between 
dislocations. In the local limit x s I, uo and I are sharply selected on a time scale of 
one rotation, and the ridge winds itself into a classic Archimedian spiral (i.e. spiral with 
a constant I) in ~ L/l rotations. In contrast, Fig. ||] shows that, in the desorpt ion-free 
limit, the transient spiral ridge evolves extremely slowly towards an Archimedian spiral via 
a progressive reduction of the step spacing away from the core. The rotation frequency of 
the center is faster than F at the beginning, as a larger step spacing allows more adatoms 
to attach to the ridge, and then slowly approaches F. In addition, Fig. || shows that the 
shape is strongly influenced by the boundaries during the initial transient. The final spacing, 
however, is independent of the system size. In order to quantify the transient dynamics, we 
calculate the surface width 
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FIG. 2. Spiral ridge at different times t' = tF, equivalent to the number of monolayers de- 
posited, for x s = oo, Fdl/D = 1.344 x 10" 5 , and L/d = 5464. 

w(t) = -< tp(x, tf- < 4j{x, t) > 2 > 1/2 (16) 

where < / >= L~ 2 j fdxdy. In steady-state, we simply have w(t — > oo) oc L/l for L ^> I. 
We plot in Fig. |3| w(i)/(L/l) as a function of Ft/ (L/l) 3 for different system sizes. The 
data collapse remarkably onto a single curve, which shows that the time to reach steady 
state scales as (L/l) 3 . This cubic law can be understood analytically by considering the 
train dynamics in the region away from the spiral core. In this region, the effect of the 
step curvature can be neglected and a one-dimensional step train is governed by the simple 
evolution equation d t l n = {F/2)(l n -\ — l n +i), where l n is the distance between the nth and 
n + 1th steps. This set of discrete equations can be transformed into a continuum equation 
for the coarse-grained surface height h(x, t) (in units of a) of the standard conserved form 



d t h = F- d x J, (17) 

where J is the surface current. Here, J = —Dd x < u > where < u >= Fl 2 /12D = 
F/12D(d x h) 2 is the the average adatom concentration between two steps, which combined 
with Eq. [T7| yields 

d ' h = F + T2^(w^)' < 18) 

This equation has a scaling solution of the form 

h(x, t) = Ft + (L/l)~h(x/(L/l),tF/(L/l) 3 ), (19) 

where the system size drops out of the resulting equation for h. Thus, the relaxation time 
depends on the third power of the system size as observed in our simulations. 
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FIG. 3. Rescaled surface width w/{L/l) as a function of rescaled time tF/(L/l) s for four 
different system sizes. The large fluctuations at early times occur during the initial winding of 
the straight ridge into a spiral shape. 



Let us now examine how these results become modified when the assumption of local 
equilibrium at the step is relaxed by letting k be finite in Eq. ^. In this case, an approximate 
expression for I can be obtained by repeating the estimate of CC with the modified boundary 
condition u(l) = v^/k = Fl/k at the edge of the circular terrace around the dislocation 
center. This yields u(0) = (1 + AD /Ik) Fl 2 /AD, where as before I = 19d Q /u(0). So for 
fast attachment (AD /Ik <C 1) we recover the previous scaling I ~ F^ 1 ^ 3 , whereas for slow 
attachment (AD /Ik ^> 1) we obtain the scaling / ~ (dok/F) 1 ^ 2 in agreement with boundary 
integral results We performed a series of simulations with a finite k by choosing A = 1 
and observed a slow cross-over towards F~ x l 2 with increasing F, which is consistent with 
this prediction. The spiral ridge was also found to relax much slower than in the local limit 
but we did not conduct a systematic finite size scaling analysis to determine the dependence 
on L. Finally, if the present analyses are extended to the case of a finite ES barrier, one 
concludes that the scaling I ~ F" 1 ' 3 remains unchanged. In this case, however, the step train 
away from the spiral core is known to be morphologically unstable [15], and this instability 
may alter spiral growth in ways that remain to be investigated. 

Several of the present predictions should be experimentally testable. If desorption is 
negligible, one should observe a dependence of the form I ~ F~ a with a between 1/3 and 1/2. 
Moreover, by examining the ridge dynamics one can obtain information about the relative 
importance of desorption and diffusion. In particular, in experiments the distance L between 
screw dislocation centers is often in the range of five to ten I If desorption is negligible, 
the cubic dependence of the relaxation time to steady-state growth on / implies that the 
spiral ridge should reach a maximum surface slope only after a few hundred monolayers are 
deposited. 

The present phase-field approach should be generally applicable to simulate a wide range 
of mesoscopic surface growth phenomena. Interesting future prospects are to include con- 
centration fluctuations to study the crossover from spiral growth to island nucleation at 
high temperature/flux, to study the effect of anisotropy of the line tension and/or the at- 
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tachment kinetics on the growth morphology, and to include an ES barrier. Finally, elastic 
effects which have recently been shown to influence spiral growth dynamics [R| could be 
incorporated by coupling the dynamics of ip and u to the strain field. 
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from supercomputer time at NERSC. We thank R. Kohn, D. Wolf, and A. Zangwill for 
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